
rm(list=ls(all=TRUE))

set.seed(123)

library(MASS)
library(matrixStats)

setwd("replication_archive/Leadership/")


### Load the data
data <- read.table("HMS_Data.raw.txt", header=T)


### Run original model

obs_Q3 <- ifelse(is.na(data$CQ3a_Clinic_or_Hospitals)==F & is.na(data$FQ3a_Clinic)==F, 1, 0)
q3_data <- subset(data, obs_Q3 == 1)

m <- glm(CQ3a_Clinic_or_Hospitals ~ FQ3a_Clinic, family=binomial(link="logit"), data=q3_data)
summary(m)

yes <- 1/(1+exp(-(c(1,1) %*% coef(m))))
no <- 1/(1+exp(-(c(1,0) %*% coef(m))))
atepoint <- yes-no
atepoint

draw <- mvrnorm(1000, coef(m), vcov(m))
yesvec <- 1/(1+exp(-(draw %*% c(1,1))))
novec <- 1/(1+exp(-(draw %*% c(1,0))))
atepointvec <- yesvec-novec
quantile(atepointvec, c(0.025, 0.975))


	
### Run the sensitivity analysis: Scenario (a)

etavec = seq(0, 0.99, length.out=100)

results <-  NULL
for(i in 1:length(etavec)){
	res <- NULL
	count <- 1
	while(count <= 500){
		# create data given eta
		q3_data$newtreat <- NA
		q3_data$newtreat[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==0] <- rbinom(length(q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==0]), 1, etavec[i])
		q3_data$newtreat[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==1] <- rbinom(length(q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==1]), 1, 1-etavec[i])
		q3_data$newtreat[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==1] <- q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==1]
		q3_data$newtreat[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==0] <- q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==0]
		if(0 %in% c(table(q3_data$newtreat, q3_data$CQ3a_Clinic_or_Hospitals))) next	

		# run model and get ATE
		m <- glm(CQ3a_Clinic_or_Hospitals ~ newtreat, family=binomial(link="logit"), data=q3_data)
		summary(m)

		yes <- 1/(1+exp(-(c(1,1) %*% coef(m))))
		no <- 1/(1+exp(-(c(1,0) %*% coef(m))))
		atepoint <- yes-no

		draw <- mvrnorm(1000, coef(m), vcov(m))
		yesvec <- 1/(1+exp(-(draw %*% c(1,1))))
		novec <- 1/(1+exp(-(draw %*% c(1,0))))
		atepointvec <- yesvec-novec

		res <- rbind(res, c(atepoint, quantile(atepointvec, 0.025), quantile(atepointvec, 0.975)))	
		count <- count+1		
	}
	results <- rbind(results, c(etavec[i], colMedians(res)))
}

colnames(results) <- c("eta", "pointest", "lowerci95", "upperci95")
results <- as.data.frame(results)


# Figure 5(a)
quartz(type="pdf", width=5, height=5, file="output/leadership_q3_sc1.pdf")
par(mar = c(4,4,0.3,0.3), mgp=c(2.5,1,0), family="CMU Serif")
plot(1:length(results$pointest), results$pointest, type="n", ylim=c(-1, 1), xlab = expression(eta), ylab="Average Treatment Effect", xaxt="n")
polygon(c(1:length(results$pointest), rev(1:length(results$pointest))), c(results$lowerci95,rev(results$upperci95)), col="grey", border=NA)
points(1:length(results$pointest), results$pointest, type="l", lwd=3)
axis(1, at=seq(1, 100, length.out=6), labels = seq(0, 1, length.out=6), las=2)
abline(h=0, col = "black", lwd=2)
lines(c(1,1), c(results$lowerci95[1], results$upperci95[1]), lwd=3)
points(1, results$pointest[1], pch=16)
dev.off()





### Run the sensitivity analysis: Scenario (b)

etavec = seq(0, 0.99, length.out=100)

results <-  NULL
for(i in 1:length(etavec)){
	res <- NULL
	count <- 1
	while(count <= 500){
		# create data given eta
		q3_data$newtreat <- NA
		q3_data$newtreat[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==1] <- rbinom(length(q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==1]), 1, etavec[i])
		q3_data$newtreat[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==0] <- rbinom(length(q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==0]), 1, 1-etavec[i])
		q3_data$newtreat[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==1] <- q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==1 & q3_data$CQ3a_Clinic_or_Hospitals==1]
		q3_data$newtreat[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==0] <- q3_data$FQ3a_Clinic[q3_data$FQ3a_Clinic==0 & q3_data$CQ3a_Clinic_or_Hospitals==0]
		if(0 %in% c(table(q3_data$newtreat, q3_data$CQ3a_Clinic_or_Hospitals))) next	

		# run model and get ATE
		m <- glm(CQ3a_Clinic_or_Hospitals ~ newtreat, family=binomial(link="logit"), data=q3_data)
		summary(m)

		yes <- 1/(1+exp(-(c(1,1) %*% coef(m))))
		no <- 1/(1+exp(-(c(1,0) %*% coef(m))))
		atepoint <- yes-no

		draw <- mvrnorm(1000, coef(m), vcov(m))
		yesvec <- 1/(1+exp(-(draw %*% c(1,1))))
		novec <- 1/(1+exp(-(draw %*% c(1,0))))
		atepointvec <- yesvec-novec

		res <- rbind(res, c(atepoint, quantile(atepointvec, 0.025), quantile(atepointvec, 0.975)))
		count <- count+1	
	}
	results <- rbind(results, c(etavec[i], colMedians(res)))
}

colnames(results) <- c("eta", "pointest", "lowerci95", "upperci95")
results <- as.data.frame(results)


# Figure 5(b)
quartz(type="pdf", width=5, height=5, file="output/leadership_q3_sc2.pdf")
par(mar = c(4,4,0.3,0.3), mgp=c(2.5,1,0), family="CMU Serif")
plot(1:length(results$pointest), results$pointest, type="n", ylim=c(-1, 1), xlab = expression(eta), ylab="Average Treatment Effect", xaxt="n")
polygon(c(1:length(results$pointest), rev(1:length(results$pointest))), c(results$lowerci95,rev(results$upperci95)), col="grey", border=NA)
points(1:length(results$pointest), results$pointest, type="l", lwd=3)
axis(1, at=seq(1, 100, length.out=6), labels = seq(0, 1, length.out=6), las=2)
abline(h=0, col = "black", lwd=2)
lines(c(1,1), c(results$lowerci95[1], results$upperci95[1]), lwd=3)
points(1, results$pointest[1], pch=16)
dev.off()






